!
! Copyright (C) 2020-2020 Quantum ESPRESSO group
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!------------------------------------------------------------------------------
PROGRAM postahc
!------------------------------------------------------------------------------
!!
!! This program
!!   - Reads the electron-phonon quantities calculated by ph.x with the
!!     electron_phonon='ahc' option.
!!   - Calculate the phonon-induced electron self-energy in the full matrix
!!     form at a given temperature.
!!
!! Input data (namelist "input") is described in Doc/INPUT_POSTAHC.
!!
!------------------------------------------------------------------------------
  USE kinds,       ONLY : DP
  USE constants,   ONLY : ry_to_kelvin, AMU_RY, RY_TO_CMM1, RYTOEV
  USE mp,          ONLY : mp_bcast, mp_sum
  USE mp_world,    ONLY : world_comm
  USE mp_global,   ONLY : mp_startup, mp_global_end
  USE mp_pools,    ONLY : npool, nproc_pool, me_pool, intra_pool_comm
  USE io_global,   ONLY : ionode_id, ionode, stdout
  USE io_files,    ONLY : tmp_dir, prefix
  USE io_global,   ONLY : meta_ionode, meta_ionode_id, qestdin
  USE environment, ONLY : environment_start, environment_end
  USE read_namelists_module, ONLY : check_namelist_read
  USE open_close_input_file, ONLY : open_input_file
  USE run_info,    ONLY : title
  USE klist,       ONLY : nks, xk
  USE symm_base,   ONLY : s, invs, nsym
  USE cell_base,   ONLY : at, bg
  USE ions_base,   ONLY : ntyp => nsp, nat, ityp
  USE wvfct,       ONLY : nbnd
  !
  IMPLICIT NONE
  !
  INTEGER, PARAMETER :: ntypx = 10
  !! Max number of atom types
  !
  ! Input variables
  !
  CHARACTER(LEN=256) :: outdir
  !! Directory containing input, output, and scratch files
  LOGICAL :: skip_upper
  !! Skip calculation of upper Fan and Debye-Waller self-energy
  LOGICAL :: skip_dw
  !! Skip calculation of Debye-Waller self-energy
  LOGICAL :: truncate_fan
  !! Truncate Fan self energy to states inside the AHC window
  LOGICAL :: truncate_dw
  !! Truncate Debye-Waller matrix element to states inside the AHC window
  LOGICAL :: adiabatic
  !! Use adiabatic approximation (i.e. ignore the phonon frequency in the denominator)
  LOGICAL :: use_irr_q
  !! Use q points from the irreducible Brillouin zone.
  INTEGER :: ahc_nbnd
  !! Number of bands for which the electron self-energy is to be computed.
  INTEGER :: ahc_nbndskip
  !! Number of bands to exclude when computing the self-energy. The
  !! self-energy is computed for ibnd from ahc_nbndskip + 1
  !! to ahc_nbndskip + ahc_nbnd.
  CHARACTER(LEN=256) :: filename
  !! Name of binary files
  CHARACTER(LEN=256) :: ahc_dir
  !! Directory where the output binary files for AHC e-ph coupling are located
  CHARACTER(LEN=256) :: flvec
  !! output file for normalized phonon displacements generated by matdyn.x
  !! The normalized phonon displacements are the eigenvectors divided by the
  !! square root of the mass then normalized. As such they are not orthogonal.
  REAL(DP) :: ahc_win_min_eV
  !! Lower bound of AHC window for the lower Fan term.
  REAL(DP) :: ahc_win_max_eV
  !! Upper bound of AHC window for the lower Fan term.
  REAL(DP) :: eta_eV
  !! Infinitesimal to add in the denominator for self-energy, in eV
  REAL(DP) :: temp_kelvin
  !! Temperature at which the electron self-energy is calculated, in Kelvins.
  REAL(DP) :: efermi_eV
  !! Fermi level energy, in eV.
  REAL(DP) :: amass_amu(ntypx)
  !! Mass of atoms one for each type, in atomic mass unit.
  !
  ! Local variables
  !
  LOGICAL :: lgamma
  !! .true. if q == Gamma
  LOGICAL :: needwf = .FALSE.
  !! do not need to read wavefunction data
  INTEGER :: ios
  !! io status
  INTEGER :: iat
  !! Counter for atoms
  INTEGER :: it
  !! Counter for types
  INTEGER :: ibnd
  !! Counter for bands
  INTEGER :: jbnd
  !! Counter for bands
  INTEGER :: ik
  !! Counter for k points
  INTEGER :: jk
  !! Counter for k points
  INTEGER :: iq
  !! Counter for q points
  INTEGER :: imode
  !! Counter for modes
  INTEGER :: rest
  !! Temporary variable for distributing q points
  INTEGER :: nq_loc
  !! Number of q points in this core
  INTEGER :: iq_from
  !! Index of first q point in this core
  INTEGER :: iq_to
  !! Index of last q point in this core
  INTEGER :: iun
  !! Unit for reading file
  INTEGER :: recl
  !! Record length for reading file
  INTEGER :: count
  !! Counter for degeneracy
  INTEGER :: nmodes
  !! Number of modes. 3 * nat
  INTEGER :: nq
  !! Number of q points
  INTEGER :: nqs
  !! Number of q points in the star
  INTEGER :: ikirr
  !! Index of irreducible k point
  INTEGER :: nkirr
  !! Number of irreducible k points
  INTEGER :: isk
  !! Index of the symmetry-unfolded k point
  INTEGER :: imq
  !! index of -q in the star of q (0 if not present)
  INTEGER :: isq (48)
  !! index of q in the star of a given sym.op
  REAL(DP) :: temperature
  !! temp_kelvin transformed from Kelvin to Ry
  REAL(DP) :: unorm
  !! Norm of u multiplied with amass
  REAL(DP) :: rval
  !! Temporary real variables
  REAL(DP) :: omega_zero_cutoff = 1.d-4
  !! Cutoff of phonon frequency. Modes with omega smaller than
  !! omega_zero_cutoff is neglected.
  REAL(DP) :: e_degen_cutoff = 2.d-5
  !! degeneracy cutoff. Ignore couping between degenerate states with energy
  !! difference less than e_degen_cutoff.
  REAL(DP) :: deltak(3)
  !! k point vector
  REAL(DP) :: sxq(3, 48)
  !! list of vectors in the star of q
  REAL(DP) :: eta
  !! eta_eV converted to Ry units
  REAL(DP) :: efermi
  !! efermi_eV converted to Ry units
  REAL(DP) :: ahc_win_min
  !! ahc_win_min_eV converted to Ry units
  REAL(DP) :: ahc_win_max
  !! ahc_win_max_eV converted to Ry units
  COMPLEX(DP) :: selfen_avg_temp(5)
  !! Diagonal self-energy averaged over degenerate states
  INTEGER, ALLOCATABLE :: ik_to_ikirr(:)
    !! Index of irreducible k point for each of the k point
  REAL(DP), ALLOCATABLE :: inv_omega(:)
  !! (nmodes) 1 / omega
  REAL(DP), ALLOCATABLE :: occph(:)
  !! (nmodes) Bose-Einstein occupation of phonon
  REAL(DP), ALLOCATABLE :: wtq(:)
  !! (nq) Weight of q points. Set to 1/nq.
  REAL(DP), ALLOCATABLE :: amass(:)
  !! Mass of atoms in Ry
  REAL(DP), ALLOCATABLE :: xq(:, :)
  !! (3, nq) q point vectors in Cartesian basis
  REAL(DP), ALLOCATABLE :: omega(:, :)
  !! (nmodes, nq) Phonon frequency
  REAL(DP), ALLOCATABLE :: etk(:, :)
  !! (nbnd, nks) Energy at k
  REAL(DP), ALLOCATABLE :: etk_all(:, :)
  !! (ahc_nbnd, nks) Energy at k
  COMPLEX(DP), ALLOCATABLE :: u(:, :, :)
  !! (nmodes, nmodes, nq) Phonon modes
  COMPLEX(DP), ALLOCATABLE :: ahc_dw(:, :, :, :, :)
  !! Debye-Waller matrix element
  COMPLEX(DP), ALLOCATABLE :: ahc_dw_trunc(:, :, :, :, :)
  !! Debye-Waller matrix element computed with truncation
  COMPLEX(DP), ALLOCATABLE :: selfen_updw(:, :, :)
  !! Debye-Waller self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_lodw(:, :, :)
  !! Debye-Waller self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_upfan(:, :, :)
  !! Upper Fan self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_lofan(:, :, :)
  !! Lower Fan self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_tot(:, :, :)
  !! Total self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_diag(:, :)
  !! Diagonal self-energy
  COMPLEX(DP), ALLOCATABLE :: selfen_diag_avg(:, :, :)
  !! Diagonal self-energy averaged over degenerate states
  !
  LOGICAL, EXTERNAL  :: imatches
  CHARACTER(LEN=256), EXTERNAL :: trimcheck
  CHARACTER(len=6), EXTERNAL :: int_to_char
  REAL(DP), EXTERNAL :: wgauss
  !
  ! ---------------------------------------------------------------------------
  !
  NAMELIST / input / prefix, outdir, skip_upper, skip_dw, ahc_nbnd, &
      ahc_nbndskip, ahc_dir, flvec, eta_eV, temp_kelvin, efermi_eV, amass_amu,  &
      ahc_win_min_eV, ahc_win_max_eV, use_irr_q, adiabatic
  !
  CALL mp_startup()
  CALL environment_start('POSTAHC')
  !
  ! ---------------------------------------------------------------------------
  !
  ! Reading input
  !
  IF (meta_ionode) THEN
    !
    ! ... Input from file (ios=0) or standard input (ios=-1) on unit "qestdin"
    !
    ios = open_input_file (  )
    !
    ! ... Read the first line of the input file
    !
    IF ( ios <= 0 ) READ( qestdin, '(A)', IOSTAT = ios ) title
    !
  ENDIF
  !
  CALL mp_bcast(ios, meta_ionode_id, world_comm )
  CALL errore( 'postahc', 'reading input file ', ABS( ios ) )
  CALL mp_bcast(title, meta_ionode_id, world_comm  )
  !
  ! Rewind the input if the title is actually the beginning of inputph namelist
  !
  IF( imatches("&input", title) ) THEN
    WRITE(stdout,'(6x,a)') "Title line not specified: using 'default'."
    title = 'default'
    IF (meta_ionode) REWIND(qestdin, iostat=ios)
    CALL mp_bcast(ios, meta_ionode_id, world_comm  )
    CALL errore('postahc', 'Title line missing from input.', abs(ios))
  ENDIF
  !
  ! Set default values for variables in namelist
  !
  prefix = 'pwscf'
  CALL get_environment_variable('ESPRESSO_TMPDIR', outdir)
  IF ( TRIM( outdir ) == ' ' ) outdir = './'
  skip_upper = .FALSE.
  skip_dw = .FALSE.
  use_irr_q = .FALSE.
  ahc_nbnd = -1
  ahc_nbndskip = 0
  ahc_dir = ' '
  flvec = ' '
  eta_eV = -9999.d0
  temp_kelvin = -9999.d0
  efermi_eV = -9999.d0
  amass_amu(:) = -9999.d0
  ahc_win_min_eV = -1.d8
  ahc_win_max_eV = 1.d8
  adiabatic = .FALSE.
  !
  ! Read the namelist input
  !
  IF (meta_ionode) THEN
     READ( qestdin, INPUT, IOSTAT = ios )
  ENDIF
  CALL check_namelist_read(ios, qestdin, "input")
  !
  CALL mp_bcast(ios, meta_ionode_id, world_comm )
  IF ( ABS(ios) /= 0 ) THEN
    CALL errore( 'postahc', 'reading input namelist', ABS( ios ) )
  ENDIF
  !
  ! Broadcast input arguments
  !
  CALL mp_bcast(prefix, ionode_id, world_comm)
  CALL mp_bcast(outdir, ionode_id, world_comm)
  CALL mp_bcast(skip_upper, ionode_id, world_comm)
  CALL mp_bcast(skip_dw, ionode_id, world_comm)
  CALL mp_bcast(use_irr_q, ionode_id, world_comm)
  CALL mp_bcast(ahc_nbnd, ionode_id, world_comm)
  CALL mp_bcast(ahc_nbndskip, ionode_id, world_comm)
  CALL mp_bcast(ahc_dir, ionode_id, world_comm)
  CALL mp_bcast(flvec, ionode_id, world_comm)
  CALL mp_bcast(eta_eV, ionode_id, world_comm)
  CALL mp_bcast(temp_kelvin, ionode_id, world_comm)
  CALL mp_bcast(efermi_eV, ionode_id, world_comm)
  CALL mp_bcast(amass_amu, ionode_id, world_comm)
  CALL mp_bcast(ahc_win_min_eV, ionode_id, world_comm)
  CALL mp_bcast(ahc_win_max_eV, ionode_id, world_comm)
  CALL mp_bcast(adiabatic, ionode_id, world_comm)
  !
  tmp_dir = trimcheck(outdir)
  ahc_dir = trimcheck(ahc_dir)
  !
  truncate_dw = skip_upper
  truncate_fan = skip_upper
  temperature = temp_kelvin / ry_to_kelvin
  ahc_win_min = ahc_win_min_eV / RYTOEV
  ahc_win_max = ahc_win_max_eV / RYTOEV
  eta = eta_eV / RYTOEV
  efermi = efermi_eV / RYTOEV
  !
  ! Read xml data. Do not need wavefunction information.
  !
  needwf = .FALSE.
  CALL read_file_new(needwf)
  !
  ! Check input argument validity
  !
  IF (ahc_nbnd < 0) CALL errore('postahc', 'ahc_nbnd must be specified', 1)
  IF (ahc_dir == '') CALL errore('postahc', 'ahc_dir must be specified', 1)
  IF (flvec == '') CALL errore('postahc', 'flvec must be specified', 1)
  IF (eta_eV < -9990.d0) CALL errore('postahc', 'eta must be specified', 1)
  IF (efermi_eV < -9990.d0) CALL errore('postahc', 'efermi must be specified', 1)
  IF (temp_kelvin < -9990.d0) CALL errore('postahc', &
      'temp_kelvin must be specified', 1)
  DO it = 1, ntyp
    IF (amass_amu(it) < -9990.d0) CALL errore('postahc', &
      'amass_amu(it) must be specified for it = 1 to ntyp', 1)
  ENDDO
  !
  ! Default parallelization is over q points. Pools not implemented.
  !
  IF (npool > 1) CALL errore('postahc', 'pools not implemented', npool)
  !
  IF (truncate_fan .neqv. truncate_dw) THEN
    WRITE(stdout, '(5x, a)') "WARNING: For double-grid calculations, it is strongly advised"
    WRITE(stdout, '(5x, a)') "to set truncate_fan and truncate_dw to the same value because"
    WRITE(stdout, '(5x, a)') "otherwise the double-grid result converge much slowly."
    WRITE(stdout, *) ""
  ENDIF
  !
  IF (adiabatic) THEN
    WRITE(stdout, '(5x, a)') "WARNING: Using the adiabatic approximation. This approximation"
    WRITE(stdout, '(5x, a)') "is inaccurate and even divergent in some materials."
    WRITE(stdout, '(5x, a)') "(See S. Poncé et al., J. Chem. Phys. 143, 102813 (2015) for details.)"
    WRITE(stdout, '(5x, a)') "One should use adiabatic = .TRUE. (default) to get accurate results."
    WRITE(stdout, *) ""
  ENDIF
  !
  ! Setup local variables, unit converstion
  !
  nmodes = 3 * nat
  !
  ALLOCATE(amass(nmodes))
  ALLOCATE(occph(nmodes))
  ALLOCATE(inv_omega(nmodes))
  ALLOCATE(etk_all(nbnd, nks))
  ALLOCATE(etk(ahc_nbnd, nks))
  ALLOCATE(ahc_dw(ahc_nbnd, ahc_nbnd, nmodes, 3, nks))
  ALLOCATE(ahc_dw_trunc(ahc_nbnd, ahc_nbnd, nmodes, 3, nks))
  ALLOCATE(selfen_updw(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(selfen_lodw(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(selfen_upfan(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(selfen_lofan(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(selfen_tot(ahc_nbnd, ahc_nbnd, nks))
  !
  etk_all = 0.d0
  etk = 0.d0
  ahc_dw = (0.d0, 0.d0)
  ahc_dw_trunc = (0.d0, 0.d0)
  selfen_updw = (0.d0, 0.d0)
  selfen_lodw = (0.d0, 0.d0)
  selfen_upfan = (0.d0, 0.d0)
  selfen_lofan = (0.d0, 0.d0)
  selfen_tot = (0.d0, 0.d0)
  !
  DO iat = 1, nat
    amass((iat-1) * 3 + 1:(iat-1) * 3 + 3) = amass_amu(ityp(iat)) * AMU_RY
  ENDDO
  !
  ! ---------------------------------------------------------------------------
  !
  ! Calculate the map of the full k points to the irreducible k points
  !
  IF (use_irr_q) CALL setup_ik_to_ikirr_mapping()
  !
  ! ---------------------------------------------------------------------------
  !
  ! Read flvec file
  !
  IF (ionode) THEN
    CALL read_flvec(flvec, nmodes, nq, xq, omega, u)
    !
    omega = omega / RY_TO_CMM1
    !
    DO iq = 1, nq
      DO imode = 1, nmodes
        unorm = SUM(CONJG(u(:, imode, iq)) * u(:, imode, iq) * amass(:))
        u(:, imode, iq) = u(:, imode, iq) / SQRT(unorm)
      ENDDO
    ENDDO
  ENDIF ! ionode
  !
  CALL mp_bcast(nq, ionode_id, world_comm)
  !
  IF (.NOT. ionode) THEN
    ALLOCATE(xq(3, nq))
    ALLOCATE(omega(nmodes, nq))
    ALLOCATE(u(nmodes, nmodes, nq))
  ENDIF
  !
  CALL mp_bcast(xq, ionode_id, world_comm)
  CALL mp_bcast(omega, ionode_id, world_comm)
  CALL mp_bcast(u, ionode_id, world_comm)
  !
  ALLOCATE(wtq(nq))
  wtq = 1.d0 / REAL(nq, KIND=DP)
  !
  ! ---------------------------------------------------------------------------
  !
  ! Distribute q points.
  !
  nq_loc = nq / nproc_pool
  rest = nq - nq_loc * nproc_pool
  IF (me_pool < rest) THEN
    nq_loc = nq_loc + 1
    iq_from = me_pool * nq_loc + 1
    iq_to = iq_from + nq_loc - 1
  ELSE
    iq_from = rest * (nq_loc + 1) + (me_pool - rest) * nq_loc + 1
    iq_to = iq_from + nq_loc - 1
  ENDIF
  WRITE(stdout, '(5x, a, I8, a)') "There are   ", nq, " q points in total."
  WRITE(stdout, '(5x, a, I8, a)') "Calculating ", nq_loc, " q points in the root core."
  !
  ! ---------------------------------------------------------------------------
  !
  ! Calculate weight of the q-points for the irreducible grid case.
  !
  IF (use_irr_q) THEN
    IF (ionode) THEN
      DO iq = 1, nq
        CALL star_q(xq(:, iq), at, bg, nsym, s, invs, nqs, sxq, isq, imq, .FALSE.)
        wtq(iq) = REAL(nqs, DP)
        IF (imq == 0) wtq(iq) = REAL(2 * nqs, DP)
      ENDDO
      WRITE(stdout, '(5x,a)') 'Check whether the sum of the weight is equal&
          & to the number of q points in the full q-point grid.'
      WRITE(stdout, '(5x,a,F15.4)') 'SUM(wtq) = ', SUM(wtq)
      wtq = wtq / SUM(wtq)
    ENDIF
    CALL mp_bcast(wtq, ionode_id, world_comm)
    !
  ENDIF ! use_irr_q
  !
  ! ---------------------------------------------------------------------------
  !
  WRITE(stdout, '(5x, a, I8)') "Total number of bands : ", nbnd
  WRITE(stdout, '(5x, "Bands to calculate the self-energy :", I8, " to ", I8)') &
    ahc_nbndskip + 1, ahc_nbndskip + ahc_nbnd
  !
  ! ---------------------------------------------------------------------------
  !
  ! Read ahc_etk_iq1 (All other q points should contain the same data.)
  !
  filename = TRIM(ahc_dir) // 'ahc_etk_iq1.bin'
  !
  INQUIRE(IOLENGTH=recl) etk_all
  OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
    ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
  READ(iun, REC=1) etk_all
  CLOSE(iun)
  !
  ! Skip ahc_nbndskip bands in etk
  !
  etk = etk_all(ahc_nbndskip+1:ahc_nbndskip+ahc_nbnd, :)
  !
  ! Read ahc_dw which does not depend on iq.
  !
  IF (ionode) THEN
    !
    CALL compute_ahc_dw_with_truncation(ahc_dw_trunc)
    !
    IF (.NOT. truncate_dw) THEN
      filename = TRIM(ahc_dir) // 'ahc_dw.bin'
      INQUIRE(IOLENGTH=recl) ahc_dw
      OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
        ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
      IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
      READ(iun, REC=1) ahc_dw
      CLOSE(iun)
      !
      ahc_dw = ahc_dw - ahc_dw_trunc
      !
    ENDIF
    !
  ENDIF ! ionode
  !
  CALL mp_bcast(ahc_dw, ionode_id, world_comm)
  CALL mp_bcast(ahc_dw_trunc, ionode_id, world_comm)
  !
  ! ---------------------------------------------------------------------------
  !
  ! Loop over q points and calculate self-energy
  !
  WRITE(stdout, *)
  WRITE(stdout,'(5x,a,I8,a)') 'Calculating electron self-energy. Loop over ', &
                              iq_to - iq_from + 1, ' q points'
  !
  DO iq = iq_from, iq_to
    !
    ! Print progress
    !
    IF (MOD(iq, 10) == 0) THEN
      WRITE(stdout, '(i8)', ADVANCE='no') iq
      IF(MOD(iq, 100) == 0 ) WRITE(stdout,*)
      FLUSH(stdout)
    ENDIF
    !
    lgamma = .FALSE.
    IF ( ALL( ABS(xq(:, iq)) < 1.d-4 ) ) lgamma = .TRUE.
    !
    ! Set Bose-Einstein occupation of phonons and set inv_omega
    !
    DO imode = 1, nmodes
      IF (omega(imode, iq) < omega_zero_cutoff) THEN
        inv_omega(imode) = 0.d0
        occph(imode) = 0.d0
      ELSE
        inv_omega(imode) = 1.d0 / omega(imode, iq)
        IF (temperature < 1.d-4) THEN
          occph = 0.d0
        ELSE
          rval = wgauss( omega(imode, iq) / temperature, -99 )
          occph(imode) = 1.d0 / (4.d0 * rval - 2.d0) - 0.5d0
        ENDIF
      ENDIF
    ENDDO
    !
    IF (.NOT. skip_dw) THEN
      CALL calc_debye_waller(iq, selfen_lodw, .TRUE.)
      IF (.NOT. truncate_dw) CALL calc_debye_waller(iq, selfen_updw, .FALSE.)
    ENDIF
    !
    IF (.NOT. truncate_fan) CALL calc_upper_fan(iq, selfen_upfan)
    !
    CALL calc_lower_fan(iq, selfen_lofan)
    !
  ENDDO ! iq
  !
  ! ---------------------------------------------------------------------------
  !
  ! Collect self-energy from all q points
  !
  CALL mp_sum(selfen_lofan, intra_pool_comm)
  CALL mp_sum(selfen_upfan, intra_pool_comm)
  CALL mp_sum(selfen_lodw, intra_pool_comm)
  CALL mp_sum(selfen_updw, intra_pool_comm)
  !
  ! If the energy is outside the AHC window, set all values to zero because
  ! the results (in particular the upper Fan self-energy) are meaningless
  DO ik = 1, nks
    DO ibnd = 1, ahc_nbnd
      IF (etk(ibnd, ik) < ahc_win_min .OR. etk(ibnd, ik) > ahc_win_max) THEN
        selfen_lofan(:, ibnd, ik) = (0.d0, 0.d0)
        selfen_lofan(ibnd, :, ik) = (0.d0, 0.d0)
        selfen_upfan(:, ibnd, ik) = 0.d0
        selfen_upfan(ibnd, :, ik) = 0.d0
        selfen_updw(:, ibnd, ik) = 0.d0
        selfen_updw(ibnd, :, ik) = 0.d0
        selfen_lodw(:, ibnd, ik) = 0.d0
        selfen_lodw(ibnd, :, ik) = 0.d0
      ENDIF
    ENDDO
  ENDDO
  !
  selfen_tot = selfen_lodw + selfen_updw + selfen_lofan + selfen_upfan
  !
  ! Write self-energy to stdout
  !
  IF (ionode) THEN
    !
    WRITE(stdout, *)
    WRITE(stdout, *)
    IF (skip_dw) THEN
      WRITE(stdout, '(5x,a)') 'skip_dw = .true.: Debye-Waller self-energy &
                              &is set to zero'
    ENDIF
    IF (skip_upper) THEN
      WRITE(stdout, '(5x,a)') 'skip_upper = .true.: the upper Fan self-energy is set &
                              &to zero, and the'
      WRITE(stdout, '(5x,a)') '                     Debye-Waller self-energy is truncated.'
    ENDIF
    !
    WRITE(stdout, *)
    WRITE(stdout, '(5x,a)') 'Diagonal electron self-energy in eV'
    WRITE(stdout, '(5x,a)') 'Self-energy of degenerate states are averaged.'
    WRITE(stdout, '(5x,a)') 'The DW and Upper_Fan terms are real-valued by construction.'
    WRITE(stdout, '(5x,a)') 'For states with energy outside the AHC window, all output are zero.'
    WRITE(stdout, '(5x,a)') 'Total = Upper_DW + Lower_DW + Upper_Fan + Lower_Fan'
    WRITE(stdout, *)
    WRITE(stdout, '(5x,a)') 'Begin postahc output'
    WRITE(stdout, '(5x,a)') '    ik  ibnd    Re[Total]     Upper_DW      Lower_DW&
                            &   Upper_Fan Re[Lower_Fan]   Im[Total]'
    !
    ALLOCATE(selfen_diag(ahc_nbnd, 5))
    ALLOCATE(selfen_diag_avg(ahc_nbnd, 5, nks))
    !
    DO ik = 1, nks
      DO ibnd = 1, ahc_nbnd
        selfen_diag(ibnd, 1) = selfen_tot(ibnd, ibnd, ik)
        selfen_diag(ibnd, 2) = selfen_updw(ibnd, ibnd, ik)
        selfen_diag(ibnd, 3) = selfen_lodw(ibnd, ibnd, ik)
        selfen_diag(ibnd, 4) = selfen_upfan(ibnd, ibnd, ik)
        selfen_diag(ibnd, 5) = selfen_lofan(ibnd, ibnd, ik)
      ENDDO
      !
      ! Average over degenerate states
      !
      DO ibnd = 1, ahc_nbnd
        !
        selfen_avg_temp = (0.d0, 0.d0)
        count = 0
        !
        DO jbnd = 1, ahc_nbnd
          !
          IF (ABS(etk(ibnd, ik) - etk(jbnd, ik)) < e_degen_cutoff) THEN
            count = count + 1
            selfen_avg_temp = selfen_avg_temp + selfen_diag(jbnd, :)
          ENDIF
          !
        ENDDO
        !
        selfen_diag_avg(ibnd, :, ik) = selfen_avg_temp / REAL(count, DP)
        !
      ENDDO ! ibnd
    ENDDO ! ik
    !
    ! Write averaged self-energy
    !
    IF (use_irr_q) THEN
      !
      ! Average over symmetry-unfolded k points
      !
      DO ikirr = 1, nkirr
        selfen_diag = (0.d0, 0.d0)
        !
        count = 0
        DO ik = 1, nks
          IF (ik_to_ikirr(ik) == ikirr) THEN
            selfen_diag = selfen_diag + selfen_diag_avg(:, :, ik)
            count = count + 1
          ENDIF
        ENDDO ! ik
        !
        selfen_diag = selfen_diag / REAL(count, KIND=DP)
        !
        DO ibnd = 1, ahc_nbnd
          WRITE(stdout, '(5x, 2I6, 6ES13.4)') ikirr, ibnd, &
              REAL(selfen_diag(ibnd, :))  * RYTOEV, &
              AIMAG(selfen_diag(ibnd, 1)) * RYTOEV
        ENDDO
        !
      ENDDO ! ikirr
      !
    ELSE ! .NOT. use_irr_q
      !
      DO ik = 1, nks
        DO ibnd = 1, ahc_nbnd
          WRITE(stdout, '(5x, 2I6, 6ES13.4)') ik, ibnd, &
              REAL(selfen_diag_avg(ibnd, :, ik))  * RYTOEV, &
              AIMAG(selfen_diag_avg(ibnd, 1, ik)) * RYTOEV
        ENDDO
      ENDDO ! ik
      !
    ENDIF ! use_irr_q
    !
    DEALLOCATE(selfen_diag)
    DEALLOCATE(selfen_diag_avg)
    !
    WRITE(stdout, '(5x,a)') 'End postahc output'
    WRITE(stdout, *)
    !
    IF (use_irr_q) THEN
      !
      WRITE(stdout, '(5x,a)') 'Using q points in the irreducible Brillouin zone.'
      WRITE(stdout, '(5x,a)') 'The off-diagonal self-energy matrix is not computed.'
      !
    ELSE ! use_irr_q
      !
      WRITE(stdout, '(5x,a)') 'Full off-diagonal complex self-energy &
                              &matrix is written in files'
      WRITE(stdout, '(5x,a)') 'selfen_real.dat and selfen_imag.dat. &
                              &These data can differ from'
      WRITE(stdout, '(5x,a)') 'the output above because the self-energy &
                              &of degenerate states are'
      WRITE(stdout, '(5x,a)') 'NOT averaged in the selfen_*.dat output.'
      !
      ! Write self-energy to selfen.dat
      !
      OPEN(NEWUNIT=iun, FILE='selfen_real.dat', FORM='FORMATTED')
      WRITE(iun, '(a)') '# Real part of full electron self-energy matrix &
                        &Re[sigma(ibnd, jbnd, ik)] in eV'
      WRITE(iun, '(a)') '#   ik  ibnd  jbnd       Total    Upper_DW    Lower_DW&
                        &   Upper_Fan   Lower_Fan'
      DO ik = 1, nks
        DO jbnd = 1, ahc_nbnd
          DO ibnd = 1, ahc_nbnd
            WRITE(iun, '(3I6, 5ES16.7)') ik, ibnd, jbnd, &
              REAL(selfen_tot(ibnd, jbnd, ik))   * RYTOEV, &
              REAL(selfen_updw(ibnd, jbnd, ik))  * RYTOEV, &
              REAL(selfen_lodw(ibnd, jbnd, ik))  * RYTOEV, &
              REAL(selfen_upfan(ibnd, jbnd, ik)) * RYTOEV, &
              REAL(selfen_lofan(ibnd, jbnd, ik)) * RYTOEV
          ENDDO
        ENDDO
      ENDDO
      CLOSE(iun)
      !
      OPEN(NEWUNIT=iun, FILE='selfen_imag.dat', FORM='FORMATTED')
      WRITE(iun, '(a)') '# Imaginary part of full electron self-energy matrix &
                        &Im[sigma(ibnd, jbnd, ik)] in eV'
      WRITE(iun, '(a)') '#   ik  ibnd  jbnd       Total    Upper_DW    Lower_DW&
                        &   Upper_Fan   Lower_Fan'
      DO ik = 1, nks
        DO jbnd = 1, ahc_nbnd
          DO ibnd = 1, ahc_nbnd
            WRITE(iun, '(3I6, 5ES16.7)') ik, ibnd, jbnd, &
              AIMAG(selfen_tot(ibnd, jbnd, ik))   * RYTOEV, &
              AIMAG(selfen_updw(ibnd, jbnd, ik))  * RYTOEV, &
              AIMAG(selfen_lodw(ibnd, jbnd, ik))  * RYTOEV, &
              AIMAG(selfen_upfan(ibnd, jbnd, ik)) * RYTOEV, &
              AIMAG(selfen_lofan(ibnd, jbnd, ik)) * RYTOEV
          ENDDO
        ENDDO
      ENDDO
      CLOSE(iun)
      !
    ENDIF ! use_irr_q
    !
  ENDIF ! ionode
  !
  ! ---------------------------------------------------------------------------
  !
  DEALLOCATE(amass)
  DEALLOCATE(wtq)
  DEALLOCATE(xq)
  DEALLOCATE(omega)
  DEALLOCATE(u)
  DEALLOCATE(occph)
  DEALLOCATE(inv_omega)
  DEALLOCATE(etk_all)
  DEALLOCATE(etk)
  DEALLOCATE(ahc_dw)
  DEALLOCATE(ahc_dw_trunc)
  DEALLOCATE(selfen_updw)
  DEALLOCATE(selfen_lodw)
  DEALLOCATE(selfen_upfan)
  DEALLOCATE(selfen_lofan)
  DEALLOCATE(selfen_tot)
  !
  ! ---------------------------------------------------------------------------
  !
  WRITE(stdout, *)
  WRITE(stdout, *)
  CALL print_clock('debye_waller')
  CALL print_clock('lower_fan')
  CALL print_clock('upper_fan')
  !
1001 CONTINUE
  !
  CALL environment_end('POSTAHC')
  CALL mp_global_end()
  !
CONTAINS
!
!------------------------------------------------------------------------------
INTEGER FUNCTION get_index_of_k_point(xk_want, xk, nks)
!------------------------------------------------------------------------------
!! Find give k vector xk_want from the array of k vectors xk.
!------------------------------------------------------------------------------
  USE kinds, ONLY : DP
  USE cell_base, ONLY : at
  !
  IMPLICIT NONE
  !
  INTEGER, INTENT(IN) :: nks
  REAL(DP), INTENT(IN) :: xk_want(3)
  REAL(DP), INTENT(IN) :: xk(3, nks)
  !
  INTEGER :: ik
  !! Index of k point
  LOGICAL :: found
  !! True if the k point in the star is found in the k point list
  !
  DO ik = 1, nks
    !
    ! Find ik such that xk(:, ik) = xk_want (mod G)
    deltak(:) = xk(:, ik) - xk_want
    !
    CALL cryst_to_cart(1, deltak, at, -1)
    !
    IF (ALL(ABS(deltak - NINT(deltak)) < 1.d-6)) THEN
      get_index_of_k_point = ik
      RETURN
    ENDIF
    !
  ENDDO
  !
  ! If the function has not returned, the k vector is not found. Return -1.
  !
  get_index_of_k_point = -1
  !
!------------------------------------------------------------------------------
END FUNCTION get_index_of_k_point
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE setup_ik_to_ikirr_mapping()
  !----------------------------------------------------------------------------
  !! When using irreducible q points, for each k point, its symmetry-equivalent
  !! k points must all be included.
  !! Here, find the list of symmetry-inequivalent (irreducible) k points and
  !! find the mapping from the full k points to the irreducible k points.
  !----------------------------------------------------------------------------
  !
  LOGICAL :: missing
  !! Set to true if any of the symmetry pairs is missing
  REAL(DP) :: xk_crys(3)
  !! k point in crystal coordinates
  !
  ALLOCATE(ik_to_ikirr(nks))
  ik_to_ikirr = -1
  ikirr = 0
  !
  missing = .FALSE.
  !
  IF (ionode) THEN
    !
    DO ik = 1, nks
      ! Skip if the irreducible k is already found.
      IF (ik_to_ikirr(ik) /= -1) CYCLE
      !
      ikirr = ikirr + 1
      ik_to_ikirr(ik) = ikirr
      !
      CALL star_q(xk(:, ik), at, bg, nsym, s, invs, nqs, sxq, isq, imq, .FALSE.)
      !
      IF (nqs == 1) CYCLE
      !
      ! Find other k points in the star
      !
      DO isk = 2, nqs
        jk = get_index_of_k_point(sxq(:, isk), xk, nks)
        !
        IF (jk == -1) THEN
          xk_crys = sxq(:, isk)
          CALL cryst_to_cart(1, xk_crys, at, -1)
          WRITE(stdout, '(5x,a,3F12.8)') 'missing k point (crystal): ', xk_crys
          missing = .TRUE.
        ELSE
          ik_to_ikirr(jk) = ikirr
        ENDIF
      ENDDO ! isk
      !
      IF (imq == 0) THEN
        DO isk = 1, nqs
          jk = get_index_of_k_point(-sxq(:, isk), xk, nks)
          !
          IF (jk == -1) THEN
            xk_crys = -sxq(:, isk)
            CALL cryst_to_cart(1, xk_crys, at, -1)
            WRITE(stdout, '(5x,a,3F12.8)') 'missing k point (crystal): ', xk_crys
            missing = .TRUE.
          ELSE
            ik_to_ikirr(jk) = ikirr
          ENDIF
        ENDDO ! isk
      ENDIF
      !
    ENDDO ! ik
    !
  ENDIF ! ionode
  !
  CALL mp_bcast(missing, ionode_id, world_comm)
  !
  IF (missing) THEN
    WRITE(stdout, *) ''
    WRITE(stdout, '(5x,a)') 'ERROR: k point in the star not found.'
    WRITE(stdout, '(5x,a)') 'To use use_irr_q = .true., the k point for the NSCF calculation list'
    WRITE(stdout, '(5x,a)') 'must contain all k points in the star once and only once.'
    WRITE(stdout, '(5x,a)') 'See above for the list of missing k points in crystal coordinates.'
    WRITE(stdout, '(5x,a)') 'Add these k points to the input of the NSCF calculation and rerun.'
    CALL errore('postahc', 'k point in the star not found.', 1)
  ENDIF
  !
  CALL mp_bcast(ik_to_ikirr, ionode_id, world_comm)
  nkirr = MAXVAL(ik_to_ikirr)
  !
  ! Print the list of irreducible k points
  !
  WRITE(stdout, *) ''
  WRITE(stdout, '(5x, a)') 'List of irreducible k points'
  WRITE(stdout, '(5x, a)') 'ik      xk(1:3) (cart. coord. in units 2pi/alat)'
  ikirr = 0
  DO ik = 1, nks
    IF (ik_to_ikirr(ik) == ik) THEN
      ikirr = ikirr + 1
      WRITE(stdout, '(I8,3F12.8)') ikirr, xk(:, ik)
    ENDIF
  ENDDO
  WRITE(stdout, *) ''
  !
END SUBROUTINE setup_ik_to_ikirr_mapping
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE compute_ahc_dw_with_truncation(ahc_dw_trunc)
  !----------------------------------------------------------------------------
  !! Compute ahc_dw_trunc using truncated sum over states.
  !! Include only the states inside the AHC window.
  !------------------------------------------------------------------------------
  !
  COMPLEX(DP), INTENT(INOUT) :: ahc_dw_trunc(ahc_nbnd, ahc_nbnd, nmodes, 3, nks)
  !! Debye-Waller matrix element
  !
  INTEGER :: ik, ib, jb, pb, idir, imode
  CHARACTER(LEN=256) :: filename
  !! Name of files to read
  COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
  !! electron-phonon matrix element
  COMPLEX(DP), ALLOCATABLE :: ahc_p(:, :, :, :)
  !! Momentum matrix elements
  !
  ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(ahc_p(nbnd, ahc_nbnd, 3, nks))
  !
  filename = TRIM(ahc_dir) // 'ahc_gkk_iq1.bin'
  CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
  !
  filename = TRIM(ahc_dir) // 'ahc_p.bin'
  CALL postahc_read_unformatted_file(filename, 1, ahc_p)
  !
  ahc_dw_trunc = (0.d0, 0.d0)
  DO ik = 1, nks
    DO pb = 1, nbnd
      IF (etk_all(pb, ik) >= ahc_win_min .AND. etk_all(pb, ik) <= ahc_win_max) THEN
        DO ib = 1, ahc_nbnd
          DO jb = 1, ahc_nbnd
            DO idir = 1, 3
              DO imode = 1, nmodes
                ahc_dw_trunc(ib, jb, imode, idir, ik) = ahc_dw_trunc(ib, jb, imode, idir, ik) + &
                  + (0.d0, 1.d0) * CONJG(ahc_gkk(pb, ib, imode, ik)) * ahc_p(pb, jb, idir, ik) &
                  - (0.d0, 1.d0) * CONJG(ahc_p(pb, ib, idir, ik)) * ahc_gkk(pb, jb, imode, ik)
              ENDDO
            ENDDO
          ENDDO
        ENDDO
      ENDIF
    ENDDO
  ENDDO
  !
  DEALLOCATE(ahc_gkk)
  DEALLOCATE(ahc_p)
  !
END SUBROUTINE compute_ahc_dw_with_truncation
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE calc_debye_waller(iq, selfen_dw, truncate)
  !----------------------------------------------------------------------------
  !!
  !! Compute Debye-Waller self-energy at iq
  !!
  !! Implements Eq.(8) of PHonon/Doc/dfpt_self_energy.pdf
  !!
  !! Here, the "operator-generalized acoustic sum rule" is used to represent
  !! Debye-Waller self-energy as a simple matrix element.
  !! See Eq.(13) of the following reference:
  !! Jae-Mo Lihm and Cheol-Hwan Park, Phys. Rev. B, 101, 121102(R) (2020).
  !!
  !----------------------------------------------------------------------------
  USE kinds,       ONLY : DP
  !
  LOGICAL, INTENT(IN) :: truncate
  !! If .TRUE., use ahc_dw_truncated, otherwise, use ahc_dw.
  INTEGER, INTENT(IN) :: iq
  COMPLEX(DP), INTENT(INOUT) :: selfen_dw(ahc_nbnd, ahc_nbnd, nks)
  !
  INTEGER :: imode, jmode, kmode
  !! Counter for modes
  INTEGER :: idir, jdir
  !! Counter for directions
  COMPLEX(DP), ALLOCATABLE :: selfen_dw_iq(:, :, :)
  !! Debye-Waller self-energy at iq
  COMPLEX(DP), ALLOCATABLE :: coeff_dw(:, :, :)
  !! Coefficients for Debye-Waller
  !
  CALL start_clock('debye_waller')
  !
  ALLOCATE(selfen_dw_iq(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(coeff_dw(nmodes, 3, nmodes))
  coeff_dw = (0.d0, 0.d0)
  selfen_dw_iq = (0.d0, 0.d0)
  !
  ! Set coeff_dw:
  ! coeff_dw(3*iat + idir, jdir, kmode)
  ! = 0.5 * ( CONJG(u(iat*3+idir, kmode)) * u(iat*3+jdir, kmode) ).real
  !   * (occph(kmode) + 0.5) * inv_omega(kmode)
  !
  DO kmode = 1, nmodes
    DO jdir = 1, 3
      DO idir = 1, 3
        DO iat = 1, nat
          imode = 3 * (iat - 1) + idir
          jmode = 3 * (iat - 1) + jdir
          coeff_dw(imode, jdir, kmode) = coeff_dw(imode, jdir, kmode) &
            + 0.5d0 * CONJG(u(imode, kmode, iq)) * u(jmode, kmode, iq) &
              * (occph(kmode) + 0.5d0) * inv_omega(kmode)
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  coeff_dw = REAL(coeff_dw, KIND=DP)
  !
  DO kmode = 1, nmodes
    DO jdir = 1, 3
      DO imode = 1, nmodes
        IF (truncate) THEN
          selfen_dw_iq = selfen_dw_iq + ahc_dw_trunc(:, :, imode, jdir, :) &
                                      * coeff_dw(imode, jdir, kmode)
        ELSE
          selfen_dw_iq = selfen_dw_iq + ahc_dw(:, :, imode, jdir, :) &
                                      * coeff_dw(imode, jdir, kmode)
        ENDIF
      ENDDO
    ENDDO
  ENDDO
  !
  selfen_dw = selfen_dw + selfen_dw_iq * wtq(iq)
  !
  DEALLOCATE(coeff_dw)
  DEALLOCATE(selfen_dw_iq)
  !
  CALL stop_clock('debye_waller')
  !
!------------------------------------------------------------------------------
END SUBROUTINE calc_debye_waller
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE calc_upper_fan(iq, selfen_upfan)
  !----------------------------------------------------------------------------
  !!
  !! Compute upper Fan (high-energy band contribution from Sternheimer
  !! equation) self-energy at iq
  !!
  !! Implements Eq.(11) of PHonon/Doc/dfpt_self_energy.pdf
  !!
  !----------------------------------------------------------------------------
  USE kinds,       ONLY : DP
  !
  INTEGER, INTENT(IN) :: iq
  COMPLEX(DP), INTENT(INOUT) :: selfen_upfan(ahc_nbnd, ahc_nbnd, nks)
  !
  CHARACTER(LEN=256) :: filename
  !! Name of files to read
  INTEGER :: iun
  !! Unit for reading file
  INTEGER :: recl
  !! Record length for reading file
  INTEGER :: ibnd, jbnd, ibq
  !! Counter for bands
  INTEGER :: imode, jmode, kmode
  !! Counter for modes
  REAL(DP) :: coeff
  !! real coefficient
  REAL(DP), ALLOCATABLE :: etq(:, :)
  !! (nbnd, nks) Energy at k+q
  COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
  !! electron-phonon matrix element
  COMPLEX(DP), ALLOCATABLE :: ahc_gkk_mode(:, :, :, :)
  !! electron-phonon matrix element in mode basis
  COMPLEX(DP), ALLOCATABLE :: selfen_upfan_iq(:, :, :)
  !! Upper Fan self-energy at iq
  COMPLEX(DP), ALLOCATABLE :: ahc_upfan(:, :, :, :, :)
  !! Upper Fan matrix element
  COMPLEX(DP), ALLOCATABLE :: ahc_upfan_mode(:, :, :, :)
  !! Upper Fan matrix element in mode basis
  !
  CALL start_clock('upper_fan')
  !
  ALLOCATE(selfen_upfan_iq(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(ahc_upfan(ahc_nbnd, ahc_nbnd, nmodes, nmodes, nks))
  ALLOCATE(ahc_upfan_mode(ahc_nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(ahc_gkk_mode(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(etq(nbnd, nks))
  selfen_upfan_iq = (0.d0, 0.d0)
  ahc_upfan_mode = (0.d0, 0.d0)
  ahc_gkk = (0.d0, 0.d0)
  ahc_gkk_mode = (0.d0, 0.d0)
  !
  ! Read files: ahc_etq, ahc_gkk, ahc_upfan_iq*.bin
  !
  filename = TRIM(ahc_dir) // 'ahc_gkk_iq' // TRIM(int_to_char(iq)) // '.bin'
  CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
  !
  filename = TRIM(ahc_dir) // 'ahc_etq_iq' // TRIM(int_to_char(iq)) // '.bin'
  !
  INQUIRE(IOLENGTH=recl) etq
  OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
    ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
  READ(iun, REC=1) etq
  CLOSE(iun)
  !
  filename = TRIM(ahc_dir) // 'ahc_upfan_iq' // TRIM(int_to_char(iq)) // '.bin'
  !
  INQUIRE(IOLENGTH=recl) ahc_upfan
  OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
    ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'Error reading ' // TRIM(filename), 1)
  READ(iun, REC=1) ahc_upfan
  CLOSE(iun)
  !
  ! rotate ahc_gkk from Cartesian to eigenmode basis
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      DO jmode = 1, nmodes
        ahc_gkk_mode(:, :, imode, ik) = ahc_gkk_mode(:, :, imode, ik) &
        + ahc_gkk(:, :, jmode, ik) * u(jmode, imode, iq)
      ENDDO
    ENDDO
  ENDDO
  !
  ! rotate ahc_upfan from Cartesian to eigenmode basis
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      DO kmode = 1, nmodes
        DO jmode = 1, nmodes
          ahc_upfan_mode(:, :, imode, ik) = ahc_upfan_mode(:, :, imode, ik) &
          + ahc_upfan(:, :, jmode, kmode, ik) &
          * CONJG(u(jmode, imode, iq)) * u(kmode, imode, iq)
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  !
  ! Add contribution of bands outside the AHC wannier window to the upfan matrix.
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      DO ibnd = 1, ahc_nbnd
        DO jbnd = 1, ahc_nbnd
          !
          DO ibq = 1, nbnd
            !
            ! Skip active states inside the ahc window
            IF (etq(ibq, ik) >= ahc_win_min .AND. etq(ibq, ik) <= ahc_win_max) CYCLE
            !
            IF (ABS(etk(ibnd, ik) - etq(ibq, ik)) > 1.d-8) THEN
              ahc_upfan_mode(ibnd, jbnd, imode, ik) &
              = ahc_upfan_mode(ibnd, jbnd, imode, ik) &
              + CONJG(ahc_gkk_mode(ibq, ibnd, imode, ik)) &
              * ahc_gkk_mode(ibq, jbnd, imode, ik) &
              / (etk(ibnd, ik) - etq(ibq, ik))
            ENDIF
            !
          ENDDO ! ibq
          !
        ENDDO ! jbnd
      ENDDO ! ibnd
    ENDDO ! imode
  ENDDO ! ik
  !
  ! Compute selfen_upfan_iq
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      !
      coeff = 0.5d0 * inv_omega(imode) * (occph(imode) + 0.5d0)
      !
      DO jbnd = 1, ahc_nbnd
        DO ibnd = 1, ahc_nbnd
          selfen_upfan_iq(ibnd, jbnd, ik) = selfen_upfan_iq(ibnd, jbnd, ik) &
            + ( ahc_upfan_mode(ibnd, jbnd, imode, ik) &
              + CONJG(ahc_upfan_mode(jbnd, ibnd, imode, ik)) ) &
            * coeff
        ENDDO
      ENDDO
      !
    ENDDO
  ENDDO
  !
  selfen_upfan = selfen_upfan + selfen_upfan_iq * wtq(iq)
  !
  DEALLOCATE(selfen_upfan_iq)
  DEALLOCATE(ahc_upfan)
  DEALLOCATE(ahc_upfan_mode)
  DEALLOCATE(ahc_gkk_mode)
  DEALLOCATE(ahc_gkk)
  DEALLOCATE(etq)
  !
  CALL stop_clock('upper_fan')
  !
!------------------------------------------------------------------------------
END SUBROUTINE calc_upper_fan
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE calc_lower_fan(iq, selfen_lofan)
  !----------------------------------------------------------------------------
  !!
  !! Compute lower Fan (low-energy band contribution) self-energy at iq
  !!
  !! Implements Eq.(10) of PHonon/Doc/dfpt_self_energy.pdf
  !!
  !----------------------------------------------------------------------------
  USE kinds,       ONLY : DP
  !
  INTEGER, INTENT(IN) :: iq
  COMPLEX(DP), INTENT(INOUT) :: selfen_lofan(ahc_nbnd, ahc_nbnd, nks)
  !
  CHARACTER(LEN=256) :: filename
  !! Name of ahc_gkk_iq*.bin file
  INTEGER :: iun
  !! Unit for reading file
  INTEGER :: recl
  !! Record length for reading file
  INTEGER :: ibq, ibk, jbk
  !! Counter for bands
  INTEGER :: imode, jmode
  !! Counter for modes
  REAL(DP) :: sign
  !! coefficients
  COMPLEX(DP) :: de1, de2
  !! coefficients
  REAL(DP), ALLOCATABLE :: etq(:, :)
  !! (nbnd, nks) Energy at k+q
  REAL(DP), ALLOCATABLE :: occq(:, :)
  !! (nbnd, nks) Fermi-Dirac occupation of energy at k+q
  COMPLEX(DP), ALLOCATABLE :: selfen_lofan_iq(:, :, :)
  !! Lower Fan self-energy at iq
  COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
  !! electron-phonon matrix element
  COMPLEX(DP), ALLOCATABLE :: ahc_gkk_mode(:, :, :, :)
  !! electron-phonon matrix element in mode basis
  COMPLEX(DP), ALLOCATABLE :: coeff(:, :, :, :)
  !! coefficient for lower Fan self-energy
  COMPLEX(DP), ALLOCATABLE :: coeff1(:, :, :, :)
  !! coefficient for lower Fan self-energy
  COMPLEX(DP), ALLOCATABLE :: coeff2(:, :, :, :)
  !! coefficient for lower Fan self-energy
  !
  CALL start_clock('lower_fan')
  !
  ALLOCATE(selfen_lofan_iq(ahc_nbnd, ahc_nbnd, nks))
  ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(ahc_gkk_mode(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(coeff(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(coeff1(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(coeff2(nbnd, ahc_nbnd, nmodes, nks))
  ALLOCATE(etq(nbnd, nks))
  ALLOCATE(occq(nbnd, nks))
  !
  selfen_lofan_iq = (0.d0, 0.d0)
  ahc_gkk = (0.d0, 0.d0)
  ahc_gkk_mode = (0.d0, 0.d0)
  coeff = (0.d0, 0.d0)
  coeff1 = (0.d0, 0.d0)
  coeff2 = (0.d0, 0.d0)
  etq = 0.d0
  occq = 0.d0
  !
  ! Read files: ahc_etq, ahc_gkk
  !
  filename = TRIM(ahc_dir) // 'ahc_gkk_iq' // TRIM(int_to_char(iq)) // '.bin'
  CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
  !
  filename = TRIM(ahc_dir) // 'ahc_etq_iq' // TRIM(int_to_char(iq)) // '.bin'
  !
  INQUIRE(IOLENGTH=recl) etq
  OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
    ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
  READ(iun, REC=1) etq
  CLOSE(iun)
  !
  ! Fermi-Dirac occupation at k+q
  !
  DO ik = 1, nks
    DO ibnd = 1, nbnd
      IF (temperature < 1.d-4) THEN
        IF (etq(ibnd, ik) < efermi) THEN
          occq(ibnd, ik) = 1.0d0
        ELSEIF (etq(ibnd, ik) > efermi) THEN
          occq(ibnd, ik) = 0.0d0
        ELSE
          occq(ibnd, ik) = 0.5d0
        ENDIF
      ELSE
        occq(ibnd, ik) = wgauss( (efermi - etq(ibnd, ik)) / temperature, -99 )
      ENDIF
    ENDDO
  ENDDO
  !
  ! rotate ahc_gkk from Cartesian to eigenmode basis
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      DO jmode = 1, nmodes
        ahc_gkk_mode(:, :, imode, ik) = ahc_gkk_mode(:, :, imode, ik) &
        + ahc_gkk(:, :, jmode, ik) * u(jmode, imode, iq)
      ENDDO
    ENDDO
  ENDDO
  !
  ! Compute coefficients
  !
  ! sign = +1 if etk(ibk, ik) > efermi
  !      = -1 otherwise
  !
  ! coeff1(ibq, ibk, imode, ik)
  ! = ( 1 - occq(ibq, ik) + occph(imode) )
  ! / ( etk(ibk, ik) - etq(ibq, ik) - omega(imode) + 1j * eta * sign )
  !
  ! coeff2(ibq, ibk, imode, ik)
  ! = ( occq(ibq, ik) + occph(imode) )
  ! / ( etk(ibk, ik) - etq(ibq, ik) + omega(imode) + 1j * eta * sign )
  !
  ! coeff = (coeff1 + coeff2) * 0.5 * inv_omega(imode)
  !
  coeff1 = 0.d0
  coeff2 = 0.d0
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      DO ibk = 1, ahc_nbnd
        !
        sign = 1.d0
        IF (etk(ibk, ik) < efermi) sign = -1.d0
        !
        DO ibq = 1, nbnd
          !
          ! Skip states outside the AHC window
          IF (etq(ibq, ik) < ahc_win_min .OR. etq(ibq, ik) > ahc_win_max) CYCLE
          !
          IF (adiabatic) THEN
            ! Adiabatic approximation: ignore omega in the denominator
            de1 = CMPLX(etk(ibk, ik) - etq(ibq, ik), eta * sign, KIND=DP)
            de2 = CMPLX(etk(ibk, ik) - etq(ibq, ik), eta * sign, KIND=DP)
          ELSE
            de1 = CMPLX(etk(ibk, ik) - etq(ibq, ik) - omega(imode, iq), eta * sign, KIND=DP)
            de2 = CMPLX(etk(ibk, ik) - etq(ibq, ik) + omega(imode, iq), eta * sign, KIND=DP)
          ENDIF
          !
          coeff1(ibq, ibk, imode, ik) = (1.d0 - occq(ibq, ik) + occph(imode)) / de1
          coeff2(ibq, ibk, imode, ik) = ( occq(ibq, ik) + occph(imode) ) / de2
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  !
  coeff = coeff1 + coeff2
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      coeff(:, :, imode, ik) = coeff(:, :, imode, ik) * 0.5d0 * inv_omega(imode)
    ENDDO
  ENDDO
  !
  ! Remove coupling with oneself
  !
  IF (lgamma) THEN
    DO ibk = 1, ahc_nbnd
      coeff(ibk + ahc_nbndskip, ibk, :, :) = (0.d0, 0.d0)
    ENDDO
  ENDIF
  !
  ! Remove coupling between degenerate states
  !
  IF (lgamma) THEN
    DO ik = 1, nks
      DO ibk = 1, ahc_nbnd
        DO ibq = 1, nbnd
          IF ( ABS(etk(ibk, ik) - etq(ibq, ik)) < e_degen_cutoff ) THEN
            coeff(ibq, ibk, :, ik) = (0.d0, 0.d0)
          ENDIF
        ENDDO
      ENDDO
    ENDDO
  ENDIF
  !
  ! Compute selfen_lofan_iq
  !
  DO ik = 1, nks
    DO imode = 1, nmodes
      !
      DO jbk = 1, ahc_nbnd
        DO ibk = 1, ahc_nbnd
          DO ibq = 1, nbnd
            !
            selfen_lofan_iq(ibk, jbk, ik) = selfen_lofan_iq(ibk, jbk, ik) &
                + CONJG(ahc_gkk_mode(ibq, ibk, imode, ik)) &
                * ahc_gkk_mode(ibq, jbk, imode, ik) &
                * ( coeff(ibq, ibk, imode, ik) &
                  + coeff(ibq, jbk, imode, ik) ) * 0.5d0
            !
          ENDDO
        ENDDO
      ENDDO
      !
    ENDDO
  ENDDO
  !
  selfen_lofan = selfen_lofan + selfen_lofan_iq * wtq(iq)
  !
  DEALLOCATE(coeff)
  DEALLOCATE(coeff1)
  DEALLOCATE(coeff2)
  DEALLOCATE(ahc_gkk)
  DEALLOCATE(ahc_gkk_mode)
  DEALLOCATE(selfen_lofan_iq)
  DEALLOCATE(etq)
  DEALLOCATE(occq)
  !
  CALL stop_clock('lower_fan')
  !
!------------------------------------------------------------------------------
END SUBROUTINE calc_lower_fan
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE read_flvec(flvec, nmodes, nq, xq, omega, u)
  !----------------------------------------------------------------------------
  !!
  !! Read flvec (.modes) file generated by matdyn.x
  !!
  !! Output
  !!   - nq : number of q points
  !!   - xq : list of q point vectors in Cartesian coordinate
  !!   - omega : phonon frequency in Ry
  !!   - u : phonon modes, normalized to satisfy u^dagger * M * u = identity
  !!         (Eq.(1) of PHonon/Doc/dfpt_self_energy.pdf)
  !!
  !----------------------------------------------------------------------------
  USE kinds,       ONLY : DP
  !
  CHARACTER(LEN=256), INTENT(IN) :: flvec
  !! Name of the modes file
  INTEGER, INTENT(IN) :: nmodes
  !! Number of modes
  INTEGER, INTENT(OUT) :: nq
  !! Number of q points
  REAL(DP), ALLOCATABLE, INTENT(INOUT) :: xq(:, :)
  !! q point vectors
  REAL(DP), ALLOCATABLE, INTENT(INOUT) :: omega(:, :)
  !! Phonon frequency
  COMPLEX(DP), ALLOCATABLE, INTENT(INOUT) :: u(:, :, :)
  !! Phonon modes
  !
  INTEGER :: i
  !! dummy variable for reading flvec
  REAL(DP) :: omega_
  !! dummy variable for reading flvec
  INTEGER :: iq
  !! Counter for q points
  INTEGER :: imode
  !! Counter for modes
  INTEGER :: iat
  !! Counter for atoms
  INTEGER :: nat
  !! number of atoms. nmodes / 3
  INTEGER :: iun
  !! Unit for reading flvec
  INTEGER :: ios
  !! io status
  !
  nat = nmodes / 3
  !
  ! Find out nq from flvec
  !
  OPEN(NEWUNIT=iun, FILE=TRIM(flvec), STATUS='OLD', FORM='FORMATTED', IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'problem reading flvec file ' // TRIM(flvec), 1)
  !
  nq = 0
  DO WHILE (.TRUE.)
    ! flvec has nmodes * (nat + 1) + 5 lines per q point
    DO i = 1, nmodes * (nat + 1) + 5
      READ(iun, '(a)', END=9000)
    ENDDO
    nq = nq + 1
  ENDDO
9000 CONTINUE
  !
  CLOSE(iun, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'problem closing flvec file ' // TRIM(flvec), 1)
  !
  ALLOCATE(xq(3, nq))
  ALLOCATE(omega(nmodes, nq))
  ALLOCATE(u(nmodes, nmodes, nq))
  !
  ! Read xq, omega, and u from flvec
  !
  OPEN(NEWUNIT=iun, FILE=TRIM(flvec), STATUS='OLD', FORM='FORMATTED', IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'problem reading flvec file ' // TRIM(flvec), 1)
  !
  DO iq = 1, nq
    READ(iun, '(a)')
    READ(iun, '(a)')
    READ(iun, '(1x,6x,3f15.8)') (xq(i, iq), i=1, 3)
    READ(iun, '(a)')
    DO imode = 1, nmodes
      READ(iun, 9010) i, omega_, omega(imode, iq)
      DO iat = 1, nat
        READ(iun, 9020) (u(i, imode, iq), i=(iat-1)*3+1, (iat-1)*3+3)
      ENDDO
    ENDDO
    READ(iun, '(a)')
  ENDDO
9010 format(5x, 6x, i5, 3x, f15.6, 8x, f15.6, 7x)
9020 format(1x, 1x, 3(f10.6, 1x, f10.6, 3x), 1x)
  !
  CLOSE(iun, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'problem closing flvec file ' // TRIM(flvec), 1)
  !
!------------------------------------------------------------------------------
END SUBROUTINE read_flvec
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
SUBROUTINE postahc_read_unformatted_file(filename, irec, array)
  !----------------------------------------------------------------------------
  !! Read a unformatted file to an array.
  !----------------------------------------------------------------------------
  !
  IMPLICIT NONE
  !
  CHARACTER(LEN=*), INTENT(IN) :: filename
  !! Name of the unformatted file to read
  INTEGER, INTENT(IN) :: irec
  !! Record index to read
  COMPLEX(DP), INTENT(INOUT) :: array(:, :, :, :)
  !! Output. Data read from the file.
  !
  INTEGER :: iun
  !! Unit for reading file
  INTEGER :: recl
  !! Record length for reading file
  INTEGER :: ios
  !! io status
  !
  INQUIRE(IOLENGTH=recl) array
  OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
    ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
  IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
  READ(iun, REC=irec) array
  CLOSE(iun, STATUS='KEEP')
END SUBROUTINE
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
END PROGRAM postahc
!------------------------------------------------------------------------------
